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ABSTRACT 


Several numerical methods used in the calculation of hydrodynamic 
shocks were investigated. Particular attention was given to the 
artificial viscosity approach of Von Neumann and Richtmyer and its 
application to the "PUFF" numerical scheme. The particle model 
approach of Ludford, Polachek and Seeger, and the method of Lax were 


also considered. 
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I. INTRODUCTION 


Frequently the numerical calculation of comoressible fluid flow is 
complicated by the presence of shock waves. The difficulties arise 
from the fact that shock waves propagate discontinuities in velocity, 
pressure, and other variables characterizing the fluid flow. Various 
approaches to this problem have been offered, each having its own 
desirable and undesirable characteristics. These approaches generally 
fall into two categories. The first category is related to the study 
of the viscosity of the fluid, whereas the second category is dependent 
on the conservation form of the hydrodynamic equations. The first 


cateqory of approaches will receive primary attention. 


IIT. NUMERICAL CALCULATION OF HYDRODYNAMIC SHOCKS 
USING AN ARTIFICIAL VISCOSITY FACTOR 


A. INTRODUCTION 

The equations describing perfect compressible fluid flow, in the 
presence of shock waves, produce solutions with discontinuities. 
Investigation of the physical situation shows that the true discontinu- 
ities, however, cannot occur due to the viscosity, or inner friction, 
of the fluid. These equations ignore the viscosity of the fluid and 
do not accurately represent the physical system. The addition of vis- 
cosity terms to this system of equations shows that the fluid behavior 
inside the shock region 1s nonlinear but continuous. The viscosity 
of the fluid is negligible outside of the shock and significant inside 
the shock. The original intent then, was to replace the shock reqion 
by a discontinuity and treat the shock as a two-sided boundary. The 
size of the jump discontinuity, or rather the boundary values, would 
be prescribed by the Rankine-Hugoniot equations. However this approach 
has several drawbacks. First, the presence of a discontinuity comoli- 
cates the use of a numerical scheme to solve the problem. Secondly, 
the shock wave is in motion, and hence, the boundary is moving. Thirdly, 
Since irreversible thermodynamic changes of state take place across a 
shock region, an increase in the specific entropy of the fluid must be 
added to the original jump conditions. And, lastly, this approach does 
not represent the physical situation in that there is no indication of 
the behavior of the fluid inside the shock region. 

Since the addition of the true viscosity terms severely complicates 


the system of differential equations, Von Neumann and Richtmyer [Ref. 7] 


have suggested the addition of a "“pseudoviscosity" term to the equa- 
tions of non-viscous flow, which will not complicate the system as much 
as the true term would. This term must, of course, conform to several 


restrictions. 


B. THE BASIC EQUATIONS 
The equations describing one-dimensional flow of a compressible 


fluid are as follows: 


Vix,t) = mat] an (2.1) 

Cat = an (2.2) 

sare + — eed (2.3) 

ae 4 (p+q) a ll (2.4) 
Vv _ aU 

a a (2:5) 


pe Fay (2.6) 


where x is the Lagrangean coordinate, X = X(x,t) is the Eulerian 
coordinate (i.e., X(x,t) gives the position, at time t, of a fluid 
element initially at position x), 0 (x) 1s the initial density, V is 
the specific volume, U is the fluid velocity, p is the static fluid 
pressive, E is, the internal energyeper unit mass, , is the ratio of 
specific heats (i.e., Gi v C.)s and q is the artificial viscosity. 
Equations (2.3), (2.4), and (2.5) are the equations of motion, of 
energy, and of continuity respectively. Equation (2.6) is the equation 


of state for a perfect gas. 


Ge THE@ARTIFIGIAlIRTSCOSIVy FACTOR 
The expression for q must satisfy the following requirements: 
1. Equations (2.3), (2.4), and (2.5) must possess solutions 
without discontinuities. 


2. The thickness of the shock must be everywhere of the same 
order as Ax (the length increment) used in the numerical 
calculation. 


3. The effect of q on eqns. (2.3) and (2.4) must be negligible 
outside the shock region. 


4. As dx+0, the solution must approach a state with a jump 
discontinuity prescribed by the Rankine-Hugoniot equations. 


Apparently, these requirements are not enough to uniquely define q. 

In any case the expression that Von Neumann and Richtmyer developed is 
2 

_ (p .cax) 

) 


q = | 3 (2.7) 
ft pt 


where c is a dimensionless constant near unity. By the use of equation 


(2.5) q can be written as 


— (cax)* au . Ee (2.8) 
y aX OX 





Von Neumann and Richtmyer have proven that q satisfies the above 
requirements for a particular case of steady-state plane shock. They 
conjecture, however, that the artificial viscosity approach would be 
equally suited to more complicated multi-dimensional flows. The problem 
they consider is the example of a one-dimensional shock wave separating 
two regions of constant state. This simulates the situation that occurs 
when a piston is pushed at a constant velocity into a long tube con- 
taining a fluid initially at rest. After the shock has traveled a 
sufficient distance from the initiating piston, it moves at a constant 


speed, s. In the absence of an artificial viscosity term, the specific 


volume, V, at some time t, is given by figure |. 


Since we are con- 


sidering steady-state solutions only, the solutions depend only on a 


linear combination of x and t given by 
w = X- St. 


Define 


Now 
U(x,t) +> U(w) = U(x - st) 


implies that equation (2.3) becomes 


ay _ dip + q) 
M dw dw 
since 
y ou. ¢ ah ow 3aU 
OW Po? OW 2. er 


yO hes apa 
W 


aX 


Similarly, equations (2.4) and (2.5) become 


e+ (p+q) a =O 


and 
ya. dU 
dw dw 
Then equations (2.11) and (2.13) aive 
ye GV. d(p+q) 
dw dw 


and equations (2.12) and (2.14) give 


dE , d[(ptq)V] Zac 
dw * dw +e ae 


(2.41) 


aU 


Po Ot 


d 
. iota) 


(Z2)) 


(2.13) 


(2.14) 


Ke. '5) 


Von Neumann and Richtmyer then integrated equations (2.13), (2.14), 


and (2.15) with respect to w giving 
MV +U = ©, (2. Gy) 


mV + ptq = Cy (2.17) 


E + (pta)¥ + 1/2 MVe = C (2.18) 


3 
as solutions of equations (2.13), (2.14), and (2.15) where C)> Co 
and C, are constants of integration. Let the initial and final values 


be denoted by: 
As woes VV, P+.» E+E. q>0 (ar 5) 
As Wo-0; VoVes P>Pe> E+Es, q>0 (2.20) 


Since V. and Ve are particular values of V, and Pas and Pe are particular 


values of p, they satisfy equation (2.17) giving 


¢ sf 
MCV. + pee OO =MMC 
fae 2a 
which implies that 
M¢(V.-V-) = (De-p.) (2.21) 
iv'f PpoPa/ - 


A similar argument using equation (2.18) yields 


Ee t+ Dp + W2MV,o = Cy 
-E, - pV, - 2 MV. = Cy 
from which it follows that 
(Ep-E,) = 1/2 MA(V,-Ve") + paVe - DeVe 


Then by equation (2.21) 


(E.-Ep) @= 1/2 hed Ce Saat ee 
fae ave ae Fee ae ee 
= 1/2 (pe-ps)(VatVe) + p.Vs - peVe 
= ali 2qhpEV Sepelan - PaVe - PeVe) > 
or 
(E,-E.) = 1/2 (p.+p,){V.-Ve) (are2\ 


Von Neumann and Richtmyer point out that equations (2.21) and (2.22) 

are independent of q, providing q-0 as wete, and in fact are the Hugoniot 
equations. Requirement (4) is then satisfied since it has just been 
shown that the Rankine-Hugoniot conditions are satisfied for the flow 
sufficiently far from the shock region. Requirement (4) may be 

examined in an alternate fashion. Let Z(x,t) be the solution to the 

set of equations describing non-viscous flow. I.e., Z(x,t) has a jump 
discontinuity, prescribed by the Rankine-Hugoniot equations, in the shock 
region. Now let w(x,t,4x) be the solution of the equations of viscous 
flow corresponding to a fixed q, or more accurately, a fixed ax. Then 


we require 


lititeoctet) = maneested) = Zig) (2.23) 
Ax +0 


By the very fashion in which q was introduced equation (2.23) is satis- 
fied. Note that q actually has the dimensions of pressure and enters 
equation (2.3) and (2.4) in the forms a and (ptq) respectively. 
Since q is continuous, £10 as q-0, and hence, all viscosity terms 
approach zero. Consequently, the system of equations describing viscous 
flow approaches the system describing nonviscous flow as our mesh size 


approaches Zero. 


The question naturally arises as to why Von Neumann and Richtmyer 
have created a process whereby decreasing the mesh size causes w(x,t,Ax) 
to approach a solution that does not represent the physical system. 

The answer is a heuristic one in that it must be possible to make q 
arbitrarily small to accomodate arbitrarily thin shock waves. It 
should also be noted that q is physically artificial and was created 
for numerical convenience only. 

To investigate the shape of the shock Von Neumann and Richtmyer 


consider solutions satisfying 


Cay or equivalently, V0. 2a) 


Equations (2.24) are normally the situation characterizing a shock 


moving to the right. Then equation (2.7) yields 
qv = (Mcax)* ce avy * (2.25) 


Now from equation (2.18) 


E+ W2Meve = Cy - (pq) 
E- 12Meve = Ca - (ptq)v - weve 
= C? - VE(ptq) - MOV], 
And then equation (2.17) gives 
2,2 _ 
E- W/2Mvo = C,- CV. (2.26) 
Then by equation (2.6) 
pV = (E)(y-1) 
eae 


MV™ + Co(y-1) - VO, (y-1) 


=a aN v" + CAV + C, (27) 
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where 


C, = -C, (y-1) (2.28) 

C,. = C., (y-1) (2.29) 
Then equation (2.17) yields 

qv = CoV - py - MAVe 


= fee 7 ll eye 
CoV - CV -C.-[ Gove + Mey?) 


py wv? - cv - c, . (2.30) 


Now for V = V. and V = Ves q = 0. Therefore since the right side of 


qv 


equation (2.30) is quadratic in V and vanishes for V = V, and V = Ves 


qV = 14] Mey? + (Co-C,)V-C 


a) 


2 5 
= me (vey, )(v-Ve) (2.31) 
Then equation (2.25) yields 
(Weax)? y= BL (v,-WV(ve¥4) 
(cax)* ( V5" - | (v,-v)(v-V,) (2.32) 


To solve equation (2.32), Von Neumann and Richtmyer proceed as follows: 








Let, 
V.+V V.-V 
7 os 3 = Cel _ A 
A a V es ? 3 Ay = y, 9 B os A, (2353) 
Then, 
dV eee We 
cox GE = ph y ay = (een 


Therefore, 


or 


giving 


Hence, 


where 


Finally, 


V.-V V.-V. 
CAX on _ pen [ lietali a ayl/2 CA , 5 iyt/e 


[41'/% [(A,-a) (AeA) I'/? 


cox SR = TYNE [(a Pa) 1'/ (2.34) 
Ne ie 
cax(q-) Se = (YAy/4 (Oe? , (2.35) 
0 Ao 
cox B= DOy'/* 1 - Beq!/2 (2.36) 


dB 
dw = [oy] '/" CAX [ (1-B°) [72 





2 51/2 dB 
w = [—r] CAX f—a-775- 
vee arcsin B , (2.378) 
_ Zeal fe 
We = [Sey CAX (2.38) 
V.-V 
A = fe le cate - ein “Sn ee 
0 Wo ; Wo 
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or 


— (2.40) 








V is obviously continuous and hence requirement (1) is satisfied. 

Von Neumann and Richtmyer state that Wo is a measure of shock thick- 
ness. Thus, if c 7s near unity, Wo is O(ax) and requirement (2) is 
Satisfied. 

Now taking the derivative of V with respect to w and setting it 


equal to zero gives 





V.-V 
dV Leet | W 
= — cos — = 0. (2201) 
dw 2 Wo Mr 
Hence 
cos = 0 ' (2.42) 
W 
0) 
iving w = elles W (n an integer), as points where V assumes 
g g 0 > > p 


its relative maxima and minima. Similarly, 





dV ee owt . W 
= { )—s sin— = 0 (2.43) 
aes ee 
implies 
ee 
sin Ww = Q : (2.44) 


0 


giving w= nnw, , (n an integer), as inflection points for V. Now 


for w = - 5 Wo» 


(2.45) 





For w = 


hoa 
= 
w 








And for w = 0 
V.+V 
i = o (2.47) 


Finally, since only solutions satisfying Ww > 0 were considered, V is 





non-oscillatory. As a result of equations (2.41) through (2.47) 
figure 2 represents our solution. 
Note that for this particular problem of steady-state plane shock, 


a . TT TT , ; OV _ : ; 
q=0 outside [-5 Wo» x Wol(i.e., the shock layer) since meee 0a this 


region. Normally, outside the shock region, q would be negligible 


in comparison with the static pressure P because of the factor (ax)° 


ov. 


a However, inside 


in equation (2.7) and a relatively small value for 


the shock layer q is comparable to p because of the abnormally large 


value of 3° encountered in that region. Hence requirement (3) is 


satisfied. Therefore, for this particular case of steady-state plane 
Shock, Von Neumann and Richtmyer have shown that their expression for 


q conforms to all the necessary restrictions. 


D. STABILITY OF THE DIFFERENTIAL EQUATIONS 

Our next concern will be the effect the introduction of an artificial 
viscosity term has on the entire system of differential equations. 
Before investigating the stability it should be noted that much of 
the computational work actually done will be omitted for the sake of 
brevity. Instead, reference will be made to the approach and ideas 


involved. 


On a given solution U(x,t), V(x,t), etc., a small perturbation 

SU, 6V, etc., is Superimposed. Then the system will be stable if the 
perturbation can be kept arbitrarily small for all t eal by initially 
choosing the perturbation at time Ue to be sufficiently small. There- 
fore consider the following variations: 

U+>+U + 6U 

V>V + 6V 

p> pt op 

> iy OG 
We then obtain our equations of first variation (i.e., higher order 


variations are considered to be negligible): 





U pteq 
ne = Ee : (2.48) 


Mt vept(y-1)8q] + Lypt(y-1)q] 2082 + y _ 











ot 
op E 
ine éV 0 (2.49) 
2 

— (cax) © aU au (cax)~ Ja] a(su) 
aq = NS ox Lael eY > 2 sop |e SU) (2.50) 

a0Ov) = tC SU) 
+) 2 Xx (2.51) 


Equations (2.48) through (2.51) are a set of simultaneous, linear 
differential equations in 6U, 6V, ép, and 6q. Their coefficients are 
composed of terms depending on the solution functions U, V, p, and q. 
Since U, V, p, and q are considered to be smooth, well-behaved functions 
of x and t, they will be considered to be constants in a small region. 
Equations (2.48) through (2.51) were combined into one equation and a 


separation of variable technique was employed. Using this technique 


possible solutions were of the form: 


6U = sue (cos kx + i sin kx) = sue Kxtat 
= sv e% (COS Kx tin KX) = ila 
(2.52) oe sper" Ce ee, ie spel karat 
éq = sqoe" (cos kx + 7 sin kx) = tg.e' 


where 6U,> 6Vo° Po» S45» k, and a are constants and k is real. Sub- 
stitution of equations (2.52) into equations (2.48) through (2.51) 
yields the following set of simultaneous homogeneous linear equations 
in bu,» 6Vo> SPo> and 6q,: 

RS = 0 (2,539 


where R is the following 4 x 4 matrix: 


20 6 ik ik 0 
0 (ey to] Ey) Lypaty-1) gor 
—_ (2.54) 
ici (cax)* aU 
\ceutggs au HH] ve at oak 
-ik 0 0 20 6 


and S is the following vector: 


Tkxtat 
6Ue 


sp e|kxtat 


S (2.55) 
ikxtat 


a ikxtat 
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Now if equation (2.53) is to hold for a nontrivial S, then 


DET(R) = O (2056) 


The characteristic equation is: 











2 3V oO 9V 21 9U 1 oV 2 3Ulau 2 

(ao,) Vot = 5 CoE (kcax) Ee - nae (kcax) axlox + k al ypt(y-1 Jaq] 
283 2} 0U a 2 9U jaU 2 0p _ 

+p PaPV + 2p.a°(kcax) Ie - 2 (keax)® 2 ie +k? 2P 20 © (2.57) 





Equation (2.57) establishes the relationship between a and k. By fixing 
k and examining the corresponding oa, we can investigate the behavior of 
the perturbation. It should be noted that the behavtor of the pertur- 
bation must be examined both jn the shock regions and tn the normal 
regions. In the shock regions all terms will be retained. In the 
normal regions, terms containing dissipative factors (7.e., terms con- 
taining ax) are dropped. Now we are interested only in the a's corres- 
ponding to very large k's. Hence we retain only the dominant terms, 

in a and ky, of equation (2.57). The dominant term in a is Po - 

The dominant term in k contains Ko and is therefore dominated by 


a > which is the dominant term in ak. Hence in the 





2p a" (kcax) 


shock regions, equation (2.57) reduces to: 


0 “anv + 20 O “(keax)* Be = 0, (2.58) 
giving 
2 iau 
a = Tebkeax)” fay, (2.59) 
o_V 


2] 


and in the normal regions 


koayp + 0 “anV =e). (2760) 
giving 
Z Koy 
a = PC (2.61) 
aa | 
0 


Thus, in the normal regions our original system of differential 
equations is stable, and in the shock regions, the system is asymptotically 
Stable. Von Neumann and Richtmyer also point out the terms in the 
equations of variation that lead to the dominant terms tn equation 


(2.57). In the shock regions: 








a(sU) ~ 39 (eu) 
<n aan Bee (2.62) 
where 
Z 
2(cAx ) - (2.63) 
Alea ° 
Vo, 
and in the normal regions: 
9*(8U) ~ ¢ 2 a°(8u) uae 
at’ 3x° 
where 
S,° : — (2.65) 
0 
0 


E. FINITE DIFFERENCE SCHEME 
To solve the system of differential equations several finite 


difference schemes could be used. The one Von Neumann and Richtmyer 
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offer is ingeniously simple. The central differences used are skill- 
fully staggered, taking advantage of the artificial nature of q. 

Let our rectangular mesh with increments Ax and At, and integers 
mon Pe 2, ..., 8 w= OF 1, 2, ... be contrived in the following 


fashion: 


UFarg = UL(i+1/2)ax, nat] 


V = Vi (detale/2)aom ndt egwete: (2.67) 


n 
i+1/2 


After rewriting equation (2.4) the finite difference scheme is as 


follows: 
; n n-1/2 n 
meee SG peep eo Fh-1/2 
2 At AX 
n-1/2 
go iyi (2.68) 
AX 
n+] n nt+1/2 ,nt1/2 
“1/2 Veet, De a ee 
. At AX 
| 2,,nt1/2 ,mtl/2,|  ntl1/2_,ntl/2 
me 2(eax)@(UME I! yor l/2y) yt l/e_yntl el 
2 4 n+] 
(Ax) (Viste + Viste) (2.70) 
n+] n n+] n 
ry Paste * Pasty + (y-1)0944/51 Visty2 ~ Viery2) 
2 At 
cyorl + yn )( n+] _ n ) 
fo see “eye? “Pi+1/2 ~ Pasi/2 0 (2.71) 


2At 


Since central differences were used the discretization error wil] 


be O(ax)* and O(at)¢. 


ae) 


F. STABILITY OF THE FINITE DIFFERENCE SCHEME 

Having shown the original system of differential equations is 
stable and having chosen a suitable finite difference scheme one must 
now show that the difference equations are in fact stable. A finite 
difference scheme is said to be unstable if the rounding error, intro- 
duced in approximating the numerical solution, grows exponentially with 
each iteration, making nonsense of our numerical data. This concept is 
analogous to the one for differential equations in that the rounding 
error corresponds to a small perturbation in the numerical solution. 

Von Neumann and Richtmyer have shown equattons (2.68) through 
(2.71) to be conditionally stable (i.e., stable only for certain com- 
binations of Ax and At). Hence certain restrictions will be placed 
on the choice of Ax and At. It should be further noted that these 
restrictions will not be the same in shock regions as in normal re- 
gions. As before much of the computational work actually done will be 
omitted for brevity and clarity. 

Outside shock regions the stability criteria of Von Neumann and 


Richtmyer is 


ee |. (2.72) 


Equation(2.72) is the usual stability criteria encountered when hydro- 
dynamic equations of the form of equation (2.64) are approximated by 
central differences. 


A similar analysis in the shock regions yields 


ea? coe (2.73) 
AX 
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From equations (2.38), (2.39), and (2.62) 


oO = 


scAX [ 


Va-Ve 





YW? cae 
Is)? ete Wo 


Equation (2.73) then becomes 


At 
—— <( 
sali | 
where 
jJo= 
Now let 
F( 


for -1/2 Wo <w < 1/2 Wo° 


Then, 


ys 


Now setting F' = 0 implies 





=~ S@C W tan yp jee" 


sel 
0 ) Mo 
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(2.74) 


(2.75) 


(2.76) 


(2.77) 


(2.78) 


: (2.79) 


since sec 6 1S never zero. Hence, we have obtained as a critical value 


5, WdeDS 4) 172 
(J+1) 


ag l/2 


= as 5 ee (2.80) 


Noting that as i> 5 from the left F(S—) +o, Since only one critical 
0 fe) 


point has been found, equation (2.80) must be the minimum value for 


F(=~) in the shock region (-1/2 Wy <W< n/2 Wo) Consequently, 
0 


equation (2.75) may be rewritten as 





2 
At Geai/fe J 
i = FT 2sc(J-1) i 


From equation (2.21) 


(p,-p. ) 

fay 1/2 

M = Tyo (2.82) 
Ve-V. 

Equations (2.22) and (2.6) give 


ee aa 
aT 


1/2 (pytp¢) (V;-Ve) 


a 


1/2 ps(V.-Ve) + Vic PelV.-Ve)> 
(2.83) 
so that 


Pet 1/2 pelWe-Ve) = [1/2 + Ladp,v, - Sif, (2,84) 
‘- lea yl i 1 eee 
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yielding 
pee Z _ Ply 4 
i Pe (V. Ve) [oop Vy Vel P; - (2.85) 


Then equations (2.10) and (2.21) give, after eliminating Pas 





c= 
or 
a 
1/2 p We 
- = Geese | EY Ae? 
Po af af f 
Finally 
1/2 
| 2 
5 + °5 cpa eek SOF (2.86) 


where Sof is the speed of sound behind the shock, relative to x. 


Now equations (2.86) and (2.81) give 


[ (ess iG 
At. __| (y+) (2.87) 


Since Ves and hence J, is generally unknown until the problem has been 
solved, the term inside the radical is replaced by its minimum with 


respect to J: 
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Let 


Then 


Setting 0' = 0 gives 


-1/2 


1/2(d-1) CC JAG] = L20-GP I - Lg 


(Gene 


-1 1 
1/2(J-1) (23-2) -f J-tor dd 


(J-1)? (OG) J] 


a2 (Gy) - 1] + 172(%1)= 0 


so that, 


But 


Teena 


yt) 


Lad 
I 
——! 


yt] 
y-l 





$ 


(2.88) 


“i 
iy] 


(2.89) 


(2.90) 


(2e91)) 


(2.92) 


Since J = 12° corresponds to an infinitely strong shock and J < | 


y-| 


implies that Q is imaginary. Since aa * 1, there are no critical 


Var 
-| 


points in the interval Pl i. Hence J = + is the point where 


y+3 


Q assumes its minimum since Q can be made arbitrarily large as J>1 


from the right. Hence 
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(2.93) 


1/2 


Finally substitution of (2.93) into (2.87) yields 


GA ae (2.94) 


Equation (2.94) is then a sufficient condition for stability of the 
difference equations inside the shock region. It should be noted that 
Von Neumann and Richtmyer have ignored boundary conditions in their 
stability analysis. Should derivative boundary conditions enter the 
problem an appropriate stability analysis will have to include a 
discussion of the difference equations used to approximate the boundary 


conditions. 
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IIT. INTEGRATION OF THE EQUATIONS OF TRUE VISCOUS FLOW 


An alternate method to the numerical calculation of hydrodynamic 
shocks is offered by Ludford, Polachek and Seeger [Ref. 5]. Their 
approach is similar to that of Von Neumann and Richtmyer in that their 
analysis examines the viscosity of the fluid. However, by replacing 
the fluid continuum by a particle model, they show that the true equa- 


tions of viscous flow can be numerically integrated. 


A. THE BASIC EQUATIONS 
The equations of one dimensional flow of a perfect viscous compres- 


sible fluid may be written 


DU 
P Dt 7s - 242, (37) 

Dp . _, ou 

Dt Px (3.2) 
iB = = of a> (324) 
TS ONT SMP = Be. (3.4) 

where 
_ 4 au 

oO a BH Ox > (3.5) 

S = C, log(p/p*). (3.6) 


T, o, S,; and u are the temperature, viscous stress, entropy per unit 
mass, and coefficient of viscosity respectively. D represents the rate 
of change of the quantity written after it, if we move with the gas 
particle (Lagrangean viewpoint). All other quantities were defined 


previously. Now p may be eliminated from eqn. (3.3) by use of eqns. 
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(3.4), (3.6) and (3.2) in the following manner: 


5M 2 op DS 2 _B (Ty OS 


OX Dt RT ca 
giving 
Oo x = Gy = (3.7) 

Since 

Be = Cy Be Clog(p/o¥)} = cyt BR- XP, (3.8) 
we have 
Then 

of - oy Gp = peek (3.10) 
finally giving 

De ~ 3x Loly-1) - vp (3.11) 


Note that eqn. (3.11) is independent of p. 

Now Ludford et al. approximate the motion of the fluid by the motion 
of a system of particles. They consider flow inside a tube of unit 
cross-section. The tube is considered to contain (N+2) particles, each 
of mass m, and each moving along a fixed straight line. Each particle 
then represents that mass m of gas which initially had the particle at 
its center (see figure 3). Hence the forces encountered within the 
original gas may be approximated by the interaction of forces between 
the particles. 


h 


If xe is the x coordinate of the ne particle, and a dot represents 


differentiation with respect to time, the viscous interaction between 


SH 


the A and (nti) particles is given by 
4 Tt) an 
C4, = FUL]. (3. 


as a result of eqn. (3.5). Now from eqn. (3.1) 


: Ey = Pr-y) (Ona, = cra 
mm kk) |) Oh me)! 
Nes Nes ~ Di n-% 
giving 
mx = a carey 5 a + (Ons, = Oe) 9 (3.14) 
Since 
m = of ed X11). (301s) 
The pressure interaction between the neh and (nt) eh particles is given 
by : ; 
X - x 
eel n 
eS Koay = &y LS Sn ay, 7 Dare (3.16) 


If the constants u and y are those of the original gas equation then 
ean. (3.2) is automatically satisfied. Hence eqns. (3.12) through 
(3.16) are sufficient to characterize the particle model. 

Before choosing a finite difference scheme, Ludford, Polachek, and 
Seeger nondimensionalize eqns. (3.12), (3.13), (3.14), and (3.16) by 


the following transformations: 


a sep i 
p = PoP , Poca Pp 
») ut 
O a VE 
—_, (Sale 
oT 
ee 
yp 
an 2 
Po 


RY 


giving in place of 


where primes denote 


A= (3) 


x! 


a 
0 


n+] 


ae 
VR0 4 


Ly (N+1) 1172 


B. FINITE DIFFERENCE EQUATIONS 


the original equations of narticle motion, 


differentiation with respect to T, and 


(3.18) 


(3.19) 


The finite difference scheme chosen to represent eqns. (3.18) is 


as follows: 


m+ | 
n+ 


2 


AT 


m+ ] 


n+] 


+ yl 


ce 


m+) 


m 
r n+] — Mtl) - (‘n Ny j 
m+, m+ m 


+h 


n+] n 
m-% 
ntl 
2 m 


(3.20) 


In the system, the quantity Path is the pressure between the nth 


and (nt1)th particles at time T = mAT. It should be noted that the 
above system cannot be solved explicitly since the superscripts of = 
are staggered. Hence eqns. (3.20) must be solved by an iterative pro- 
cedure at every AT/2 time increment. 

The stability analysis of the finite difference equations is quite 
similar to that of Von Neumann and Richtmyer. Equations (3.20) are 


numerically stable if, in the normal regions, 


AT < ody (3.21) 
The system is unconditionally stable in shock regions. 

The approach taken by Ludford, Polachek, and Seeger has some 
advantages over that of Von Neumann and Richtmyer. First, in the 
particle model approach, the stability requirements are somewhat less 
stringent. Secondly, this approach gives a fairly accurate representation 
of the behavior of the fluid in the shock regions, whereas the artificial 
viscosity approach does not represent the physical system in the shock 
layers. However the approach by Ludford et al. does have a major dis- 
advantage. If the fluid under investigations has a very low viscosity 
an inconveniently fine mesh will have to be used in order that the finite 
difference approximations are accurate beyond the shock front. Otherwise 
AX will be larger than the thickness of the shock wave. It should be 
recalled that in the artificial viscosity approach q automatically 


adjusted to the shock thickness. 
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IV. THE METHOD OF LAX 


It was stated earlier that approaches to the numerical calculation 
of compressible fluid flow generally fall into two categories. The 
first cateaory consisted of methods which examined the viscosity of the 
fluid. The motivation for this approach was that when a nonlinear term, 
multiplied by a small coefficient, is introduced irto a differential 
equation, it may produce large changes in the behavior of the solution. 
The second category consists of methods which examine the basic equations 
of nonviscous flow, while allowing discontinuous solutions. An ingenious 
method developed by Peter Lax [Ref. 4] belongs to the second category. 

It is evident that the approaches discussed so far are heuristic in 
nature and lack theoretical preliminaries. Lax succeeded in making 
several theoretical observations which could tie together the various 
methods discussed in this paper. Unfortunately several of his observa- 
tions have been proven only for particular cases. 

Observing that all nonlinear systems of fluid dynamics satisfy 
certain "conservation laws’, Lax states that any hydrodynamic system can 
generally be brought into the form 


ve + 8 = QO, (4.1) 
where U is a column vector of unknown functions, F is a column vector 
such that F = F(x,t,U), and B is a vector coefficient. U is said to be 
a weak solution of eqn. (4.1) with initial value © if the integral 
relation 


Oo Oo oo 


/ es {WU cael WR} dxdt + f W(x,0)O(x)dx = O (4.2) 


0 -co ~co 
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holds for every test vector W which has continuous first derivatives and 
vanishes outside of some bounded region in the x,t-plane. Equation (4.2) 


is obtained by multiplying ean. (4.1) by W and then integrating by parts 


as follows: 
Pa {WU, + WFO + WB} dxdt 
0  - ‘ 
= f [WU] dx-/s f UW dxdt SE | ait 
00 0 frac 0 -0o 
-f f Fl dxdt + f ff WBdxdt (4.3) 
0 -0O 0 =00 
ieee tema 
Now I, = 0 since W vanishes outside some bounded region. Similarly 


3 


uA Oe (4.4) 


Hence we obtain ean. (4.2). 

Clearly, the only requirement for U to be a weak solution is that 
it be continuous almost everywhere so that the (Riemann) integrals in 
ean. (5.2) exist. Hence weak solutions need not be differentiable. 
The motivation for developing the notion of weak solutions is that in 
physical systems we are concerned with discontinuous functions that 
satisfy eqn. (4.1) almost everywhere. 

The Rankine-Hugoniot equations now have an interesting interpretation 
in terms of weak solutions; if U, and Us are two genuine solutions of 
ean. (4.1) whose domains in the x,t-plane are separated by a smooth curve, 


the two taken together will constitute a weak solution if and only if the 
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slope m of the separating curve and the value of U, and Uy on the curve 


satisfy 


(U, - U ee (4.5) 


m (Uy = Up) 2) 
Equation (4.5) corresponds precisely to eqns. (2.21) and (2.22). 
Having generalized the concept of a solution to a differential 
equation, we might expect some of the properties of the original concept 
to be lost. For instance, initial values do not in general determine a 
unique weak solution. This property raises immediate difficulties. In 
nature nonviscous flow is described by a unique weak solution to the 
hvdrodynamic equations, given an initial vector. Hence if our mathemati— 
cal model is a meaninaful one, there must be some other principle that 
uniquely defines a weak solution. Lax offers some possibilities, the 
most likely being that the weak solutions occuring in nature are limits 
of viscous flows. It should be possible then to determine some relation-- 
shin between weak solutions and solutions to the viscous flow problem. 
THEOREM: 


Consider the nonlinear parabolic system 


U, + i” +B = MU (4.6) 


with initial vector U(x,0) = (x). Here AUxx corresponds to a viscosity 
factor. Given that the strong limit (as A+0) of the net of solution 


functions U,(x,t) exists and is equal to U(x,t), then U(x,t) is a weak 


, | 
solution of eqn. (4.1). 
PROOF: 

Consider an arbitrary twice differentiable test vector Wl. Multiplying 


eqn. (4.6) by ' and intearating by parts yields 


3/ 


Le @) 


; oil {WU, + WF(U, ) - WB} dxdt + JS W(x,0)6(x)dx 


=O 


ie ere) (do) 
= - Af f WU, dxdt 
ea 


Now keepina ® and W fixed, and letting A=0, the left side of eqn. (4.7) 
approaches the left side of eqn. (4.2), and the right side of eqn. (4.7) 
approaches zero. Hence U(x,t) satisfies eqn. (4.2) for all twice 
differentiable test functions. Note that U had to be a strong limit 
of U, le Ciel fi |, - U| + 0 over any bounded region in the x,t-plane) 
in order that F(U, )>F(U). If U,>U only in the weak sense “se 
U, (x,t) U(x,t) vy (x,t))there is no quarantee that F(U, ) would converge 
to F(U). Now it must be shown that eqn. (4.2) is satisfied for any 
arbitrary test vector W which is continuously once differentiable. Lax 
states that since U satisfies eqn. (4.2) for all twice differentiable 
test vectors |, a fortiore, U satisfies ean. (4.2) for all once different- 
jable test ee This author disagrees with this argument since the 
class of once differentiable functions is laraer than the class of twice 
differentiable functions. Though possibly true, the statement requires 
proof. Several approaches were attempted and proved to be unsuccessful. 
Nne apparent way to eliminate the difficulty is to restrict our attention 
to just twice differentiable test vectors in our original definition of a 
weak solution. 

Having shown for a particular case, that the limit of viscous flow 


(as \>0) is a weak solution, Lax conjectures that all viscosity methods 


Vax, P. D., "Weak Solutions of Nonlinear Hyperbolic Equations and 


Their Numerical Computation," Communications on Pure and Applied 
Mathematics, v. 7, p. 163, 1954. 


38 


Should converge to the same weak solution. However he proposes a 
difterent limiting process, namely a special finite difference scheme, 
which is independent of the viscosity of the fluid. This approach has 

a sianificant advantage over the viscosity approach in that it accurately 
represents the actual behavior of the fluid in every region of the flow. 
Lax's method replaces the differentiations by finite-difference operations 


accordina to the scheme 


] 


FL (x,t)> Tay LF(xtdx,t) ~ f(x-Ax,t) | 
(4.8) 
Fitter Le(xsteat) - Tetent) = fixcdx tds 


Finally, substantial numerical evidence supports Lax's conjecture 
that when the above finite difference scheme is applied to any single 


homogeneous first order conservation law 


u, + rane = 0 lu 0 (4.9) 


with U(x,0) = (x), the solution U(x,t,4x) approaches the same weak 


solution aenerated by the viscosity methods. The solution is given by 





ee) (4.10) 
where y= y (xst) maximizes 
y x= 
f o(s)ds + t(=) (4.11) 
0 
and 
flgicuiee= 5° 5 G(s) = g (4.12) 


In most physical situations U(x,t) is a piecewise differentiable function. 


Then it is an easy matter to show U(x,t) is a weak solution to ean. (4.9) 


Be 


Since verification of eqn. (4.2) may be avoided. All that must be shown 
is that U(x,t) satisfies ean. (4.9), wherever it has well-defined first 
derivatives. Since the first derivatives of U(x,t) are undefined at most 
on a set of measure zero, and U(x,t) satisfies eqn. (4.9) everywhere 
else, U(x,t) must by necessity satisfy the intearal relation (4.2) 
evervwhere. Hence (4.10) is a weak solution of ean. (4.9). The intearal 


eqn. (4.11) is necessary to uniquely define the weak solution. 
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V. CONCLUSIONS 


The original intention of this paper was to analyze difficulties in 
a numerical system called "PUFF" [Ref. 2]. "PUFF" is an attempt to 
numerically compute the reaction of a multi-layered medium to violent 
Shocks. The Puff code employs Von Neumann and Richtmyer's artificial 
viscosity term. It has been shown that there were very large discrep- 
ancies between the "PUFF" solution and the classical solution inside 
the shock regions [Ref. 1]. The reason for these errors should now be 
apparent. Von Neumann and Richtmyer's artificial viscosity term had 
to comply with certain requirements. However these requirements were 
to a large degree independent of actual physical considerations and, 
hence, were not sufficiently restrictive. Consequently, in shock 
regions where viscosity is significant, the artificial term may not 
accurately describe the actual fluid flow. This difficulty becomes 
critical in a multi-layered medium since the number of shock waves is 


increased when the original shock reflects from several boundaries. 
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